<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of gssm_n1</title>
  <meta name="keywords" content="gssm_n1">
  <meta name="description" content="GSSM_N1  Generalized state space model for simple nonlinear system">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../../menu.html">Home</a> &gt;  <a href="#">ReBEL-0.2.7</a> &gt; <a href="#">examples</a> &gt; <a href="#">gssm</a> &gt; gssm_n1.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../../menu.html"><img alt="<" border="0" src="../../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="menu.html">Index for .\ReBEL-0.2.7\examples\gssm&nbsp;<img alt=">" border="0" src="../../../right.png"></a></td></tr></table>-->

<h1>gssm_n1
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="box"><strong>GSSM_N1  Generalized state space model for simple nonlinear system</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="box"><strong>function [varargout] = model_interface(func, varargin) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre class="comment"> GSSM_N1  Generalized state space model for simple nonlinear system

 The model is a simple scalar nonlinear system with Gamma process and Gaussian observation noise.
   Copyright (c) Oregon Health &amp; Science University (2006)

   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for
   academic use only (see included license file) and can be obtained from
   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the
   software should contact rebel@csee.ogi.edu for commercial licensing information.

   See LICENSE (which should be part of the main toolkit distribution) for more
   detail.</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../../matlabicon.gif)">
<li><a href="../../.././ReBEL-0.2.7/core/consistent.html" class="code" title="function errstring = consistent(ds, type)">consistent</a>	CONSISTENT   Check ReBEL data structures for consistentency.</li><li><a href="../../.././ReBEL-0.2.7/core/gennoiseds.html" class="code" title="function NoiseDS = gennoiseds(ArgDS)">gennoiseds</a>	GENNOISEDS    Generates a NoiseDS data structure describing a noise source.</li></ul>
This function is called by:
<ul style="list-style-image:url(../../../matlabicon.gif)">
<li><a href="../../.././ReBEL-0.2.7/examples/state_estimation/demse2.html" class="code" title="">demse2</a>	DEMSE2  Demonstrate state estimation on a simple scalar nonlinear (time variant) problem</li></ul>
<!-- crossreference -->

<h2><a name="_subfunctions"></a>SUBFUNCTIONS <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<ul style="list-style-image:url(../../../matlabicon.gif)">
<li><a href="#_sub1" class="code">function model = init(init_args)</a></li><li><a href="#_sub2" class="code">function model = setparams(model, params, index_vector)</a></li><li><a href="#_sub3" class="code">function new_state = ffun(model, state, V, U1)</a></li><li><a href="#_sub4" class="code">function observ = hfun(model, state, N, U2)</a></li><li><a href="#_sub5" class="code">function tranprior = prior(model, nextstate, state, U1, pNoiseDS)</a></li><li><a href="#_sub6" class="code">function llh = likelihood(model, obs, state, U2, oNoiseDS)</a></li><li><a href="#_sub7" class="code">function innov = innovation(model, obs, observ)</a></li><li><a href="#_sub8" class="code">function out = linearize(model, state, V, N, U1, U2, term, index_vector)</a></li></ul>
<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre>0001 <span class="comment">% GSSM_N1  Generalized state space model for simple nonlinear system</span>
0002 <span class="comment">%</span>
0003 <span class="comment">% The model is a simple scalar nonlinear system with Gamma process and Gaussian observation noise.</span>
0004 <span class="comment">%   Copyright (c) Oregon Health &amp; Science University (2006)</span>
0005 <span class="comment">%</span>
0006 <span class="comment">%   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for</span>
0007 <span class="comment">%   academic use only (see included license file) and can be obtained from</span>
0008 <span class="comment">%   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the</span>
0009 <span class="comment">%   software should contact rebel@csee.ogi.edu for commercial licensing information.</span>
0010 <span class="comment">%</span>
0011 <span class="comment">%   See LICENSE (which should be part of the main toolkit distribution) for more</span>
0012 <span class="comment">%   detail.</span>
0013 
0014 <span class="comment">%=============================================================================================</span>
0015 
0016 <a name="_sub0" href="#_subfunctions" class="code">function [varargout] = model_interface(func, varargin)</a>
0017 
0018   <span class="keyword">switch</span> func
0019 
0020     <span class="comment">%--- Initialize GSSM data structure --------------------------------------------------------</span>
0021     <span class="keyword">case</span> <span class="string">'init'</span>
0022       model = <a href="#_sub1" class="code" title="subfunction model = init(init_args)">init</a>(varargin);
0023         error(<a href="../../.././ReBEL-0.2.7/core/consistent.html" class="code" title="function errstring = consistent(ds, type)">consistent</a>(model,<span class="string">'gssm'</span>));               <span class="comment">% check consistentency of initialized model</span>
0024       varargout{1} = model;
0025 
0026     <span class="comment">%--------------------------------------------------------------------------------------------</span>
0027     <span class="keyword">otherwise</span>
0028 
0029       error([<span class="string">'Function '''</span> func <span class="string">''' not supported.'</span>]);
0030 
0031   <span class="keyword">end</span>
0032 
0033 
0034 <span class="comment">%===============================================================================================</span>
0035 <a name="_sub1" href="#_subfunctions" class="code">function model = init(init_args)</a>
0036 
0037   model.type = <span class="string">'gssm'</span>;                     <span class="comment">% object type = generalized state space model</span>
0038   model.tag  = <span class="string">'GSSM_N1'</span>;                  <span class="comment">% ID tag</span>
0039 
0040   model.ffun      = @<a href="#_sub3" class="code" title="subfunction new_state = ffun(model, state, V, U1)">ffun</a>;                   <span class="comment">% file handle to FFUN</span>
0041   model.hfun      = @<a href="#_sub4" class="code" title="subfunction observ = hfun(model, state, N, U2)">hfun</a>;                   <span class="comment">% file handle to HFUN</span>
0042   model.prior     = @<a href="#_sub5" class="code" title="subfunction tranprior = prior(model, nextstate, state, U1, pNoiseDS)">prior</a>;
0043   model.likelihood = @<a href="#_sub6" class="code" title="subfunction llh = likelihood(model, obs, state, U2, oNoiseDS)">likelihood</a>;            <span class="comment">% file handle to LIKELIHOOD</span>
0044   model.innovation = @<a href="#_sub7" class="code" title="subfunction innov = innovation(model, obs, observ)">innovation</a>;            <span class="comment">% file handle to INNOVATION</span>
0045   model.linearize  = @<a href="#_sub8" class="code" title="subfunction out = linearize(model, state, V, N, U1, U2, term, index_vector)">linearize</a>;              <span class="comment">% file handle to LINEARIZE</span>
0046   model.setparams  = @<a href="#_sub2" class="code" title="subfunction model = setparams(model, params, index_vector)">setparams</a>;              <span class="comment">% file handle to SETPARAMS</span>
0047 
0048   model.statedim   = 1;                      <span class="comment">%   state dimension</span>
0049   model.obsdim     = 1;                      <span class="comment">%   observation dimension</span>
0050   model.paramdim   = 2;                      <span class="comment">%   parameter dimension</span>
0051   model.U1dim      = 1;                      <span class="comment">%   exogenous control input 1 dimension</span>
0052   model.U2dim      = 1;                      <span class="comment">%   exogenous control input 2 dimension</span>
0053   model.Vdim       = 1;                      <span class="comment">%   process noise dimension</span>
0054   model.Ndim       = 1;                      <span class="comment">%   observation noise dimension</span>
0055 
0056   Arg.type = <span class="string">'gamma'</span>;
0057   Arg.dim = model.Vdim;
0058   Arg.alpha = 3;
0059   Arg.beta  = 0.5;
0060   model.pNoise = <a href="../../.././ReBEL-0.2.7/core/gennoiseds.html" class="code" title="function NoiseDS = gennoiseds(ArgDS)">gennoiseds</a>(Arg);   <span class="comment">% process noise : Gamma(3,0.5) noise source</span>
0061 
0062   Arg.type = <span class="string">'gaussian'</span>;
0063   Arg.cov_type = <span class="string">'full'</span>;
0064   Arg.dim = model.Ndim;
0065   Arg.mu = 0;
0066   Arg.cov  = 1e-5;
0067   Arg.cov_type = <span class="string">'full'</span>;
0068   model.oNoise = <a href="../../.././ReBEL-0.2.7/core/gennoiseds.html" class="code" title="function NoiseDS = gennoiseds(ArgDS)">gennoiseds</a>(Arg);     <span class="comment">% observation noise : zero mean white Gaussian noise, cov=0.2</span>
0069 
0070   model.params = zeros(model.paramdim,1);
0071 
0072   model = <a href="#_sub2" class="code" title="subfunction model = setparams(model, params, index_vector)">setparams</a>(model,[4e-2 0.5]);   <span class="comment">% [omega phi]</span>
0073 
0074 
0075 <span class="comment">%===============================================================================================</span>
0076 <a name="_sub2" href="#_subfunctions" class="code">function model = setparams(model, params, index_vector)</a>
0077 
0078   <span class="keyword">if</span> (nargin==2)
0079     model.params = params(:);
0080   <span class="keyword">elseif</span> (nargin==3)
0081     model.params(index_vector) = params(:);
0082   <span class="keyword">else</span>
0083     error(<span class="string">'[ setparams ] Incorrect number of input arguments.'</span>);
0084   <span class="keyword">end</span>
0085 
0086 <span class="comment">%===============================================================================================</span>
0087 <a name="_sub3" href="#_subfunctions" class="code">function new_state = ffun(model, state, V, U1)</a>
0088 
0089   new_state      = 1 + sin(model.params(1)*pi.*U1) + model.params(2)*state;
0090 
0091   <span class="keyword">if</span> ~isempty(V)
0092     new_state = new_state + V;
0093   <span class="keyword">end</span>
0094 
0095 <span class="comment">%===============================================================================================</span>
0096 <a name="_sub4" href="#_subfunctions" class="code">function observ = hfun(model, state, N, U2)</a>
0097 
0098   [dim,nop] = size(state);
0099 
0100   observ = zeros(model.obsdim,nop);
0101 
0102   <span class="keyword">for</span> k=1:nop,
0103     <span class="keyword">if</span> (U2(k) &lt;= 30),
0104        observ(k) = model.params(2)*state(:,k).^2;
0105     <span class="keyword">else</span>
0106        observ = model.params(2)*state - 2;
0107     <span class="keyword">end</span>
0108   <span class="keyword">end</span>
0109 
0110   <span class="keyword">if</span> ~isempty(N)
0111     observ = observ + N;
0112   <span class="keyword">end</span>
0113 
0114 
0115 <span class="comment">%===============================================================================================</span>
0116 <a name="_sub5" href="#_subfunctions" class="code">function tranprior = prior(model, nextstate, state, U1, pNoiseDS)</a>
0117 
0118   X = nextstate - <a href="#_sub3" class="code" title="subfunction new_state = ffun(model, state, V, U1)">ffun</a>(model, state, [], U1);
0119 
0120   tranprior = pNoiseDS.likelihood( pNoiseDS, X);
0121 
0122 
0123 <span class="comment">%===============================================================================================</span>
0124 <a name="_sub6" href="#_subfunctions" class="code">function llh = likelihood(model, obs, state, U2, oNoiseDS)</a>
0125 
0126   X = obs - <a href="#_sub4" class="code" title="subfunction observ = hfun(model, state, N, U2)">hfun</a>(model, state, [], U2);
0127 
0128   llh = oNoiseDS.likelihood( oNoiseDS, X);
0129 
0130 
0131 
0132 <span class="comment">%===============================================================================================</span>
0133 <a name="_sub7" href="#_subfunctions" class="code">function innov = innovation(model, obs, observ)</a>
0134 
0135   innov = obs - observ;
0136 
0137 
0138 <span class="comment">%===============================================================================================</span>
0139 <a name="_sub8" href="#_subfunctions" class="code">function out = linearize(model, state, V, N, U1, U2, term, index_vector)</a>
0140 
0141   <span class="keyword">if</span> (nargin&lt;7)
0142     error(<span class="string">'[ linearize ] Not enough input arguments!'</span>);
0143   <span class="keyword">end</span>
0144 
0145   <span class="comment">%--------------------------------------------------------------------------------------</span>
0146   <span class="keyword">switch</span> (term)
0147 
0148     <span class="keyword">case</span> <span class="string">'A'</span>
0149       <span class="comment">%%%========================================================</span>
0150       <span class="comment">%%%             Calculate A = df/dstate</span>
0151       <span class="comment">%%%========================================================</span>
0152       out = model.params(2);
0153 
0154     <span class="keyword">case</span> <span class="string">'B'</span>
0155       <span class="comment">%%%========================================================</span>
0156       <span class="comment">%%%             Calculate B = df/dU1</span>
0157       <span class="comment">%%%========================================================</span>
0158       out = [];
0159 
0160     <span class="keyword">case</span> <span class="string">'C'</span>
0161       <span class="comment">%%%========================================================</span>
0162       <span class="comment">%%%             Calculate C = dh/dx</span>
0163       <span class="comment">%%%========================================================</span>
0164       <span class="keyword">if</span> (U2 &lt;= 30)
0165         out = 2*model.params(2)*state;
0166       <span class="keyword">else</span>
0167         out = model.params(2);
0168       <span class="keyword">end</span>
0169 
0170     <span class="keyword">case</span> <span class="string">'D'</span>
0171       <span class="comment">%%%========================================================</span>
0172       <span class="comment">%%%             Calculate D = dh/dU2</span>
0173       <span class="comment">%%%========================================================</span>
0174       out = [];
0175 
0176     <span class="keyword">case</span> <span class="string">'G'</span>
0177       <span class="comment">%%%========================================================</span>
0178       <span class="comment">%%%             Calculate G = df/dv</span>
0179       <span class="comment">%%%========================================================</span>
0180       out = 1;
0181 
0182     <span class="keyword">case</span> <span class="string">'H'</span>
0183       <span class="comment">%%%========================================================</span>
0184       <span class="comment">%%%             Calculate H = dh/dn</span>
0185       <span class="comment">%%%========================================================</span>
0186       out = 1;
0187 
0188     <span class="keyword">case</span> <span class="string">'JFW'</span>
0189       <span class="comment">%%%========================================================</span>
0190       <span class="comment">%%%             Calculate  = dffun/dparameters</span>
0191       <span class="comment">%%%========================================================</span>
0192       out = [cos(model.params(1)*pi*U1)*pi*U1 state];
0193 
0194 
0195     <span class="keyword">case</span> <span class="string">'JHW'</span>
0196       <span class="comment">%%%========================================================</span>
0197       <span class="comment">%%%             Calculate  = dhfun/dparameters</span>
0198       <span class="comment">%%%========================================================</span>
0199       <span class="keyword">if</span> (U2 &lt;= 30)
0200         out = [0 state^2];
0201       <span class="keyword">else</span>
0202         out = [0 state];
0203       <span class="keyword">end</span>
0204 
0205     <span class="keyword">otherwise</span>
0206       error(<span class="string">'[ linearize ] Invalid model term requested!'</span>);
0207 
0208   <span class="keyword">end</span>
0209 
0210   <span class="keyword">if</span> (nargin==8), out = out(:,index_vector); <span class="keyword">end</span>
0211 
0212   <span class="comment">%--------------------------------------------------------------------------------------</span></pre></div>
<hr><address>Generated on Tue 26-Sep-2006 10:36:21 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>